/***************************************************************************
 *   Copyright (C) 2008 by Ivan Blagoev Topolsky   *
 *   ivan.topolsky@genebio.com   *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 3 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#ifdef _MSC_VER

/**
If you are using this file from Microsoft Visual C++ (MSVC++) :

!!! vital informations are in ReadMe.txt !!!

basic project manipulation :

a MSVC's solution is equivalent to a whole package in the automake/autoconf world
a MSVC's project is roughly equivalent to a target inside a automake/autoconf project

bol-wiff: is the VC project creating the final EXE file
captainbol-lib: is the VC project creating a Microsoft .LIB with the objects from Captain-Bol
libprotobuf: is found in the vsproject subfolder of Google's protocole buffers


to add a new include path :
	- right click on the project "bol-wiff" and select "properties"
	- switch drop box "Configuration:" to "All"
	- in tree pick "Configuration Properties" -> "C/C++" -> "General"
	- paths are in "Additionnal include directories"

to add a new library path :
	- switch drop box "Configuration:" to "All"
	- in tree pick "Configuration Properties" -> "Linker" -> "General"
	- paths are in "Additionnal library directories"
*/

#endif


/*
 * precompiled stuff
 */

/// The actual include go in file STDAFX.H due to the way MSVC recompiled header works
// NOTE when adding an include somewhere, add it to "stdafx.h" otherwise the precompiled header mecanism breaks appart
#include "stdafx.h"
#pragma hdrstop
// NOTE non-precompiled things go bellow


/**
 * macros to help handling COM stuff
 */

#ifdef _MSC_VER
/**
Microsoft's COM interface :

in MSVC an '#import' directive will automatically generate and include the necessary files
(tlh / tli - header and implementation) out of the tlb (typelibrary) included in the COM .DLL

These systems declare pointer to interfaces : 
- I{class_name}Ptr which are derivative of IDisplayPtr (itself driveting from IUnkown)
- These are smart pointers
- They have a CreateInstance method which populate the pointer with a NEW instance of a class (given by GUID)
  ( be careful to type the correct "{class_name}", otherwise, one may get : 80040154 "Class not registered")
- They point to a structure with all the necessary method (inline) to call the COM's method (stored in the object's virtual method table (vptr) - the only thing portable accross objects in the Windows world)
- These methods are defined at compile time by the TLH (= early binding)
- The .DLL is required for compiling (or at least, one should include prebuilt TLH/TLI)
- This system is highly sensitive to versions (if position inside the vptr changes, everything breaks)

*/



// macros used to populate a pointer with an instance of a class
// or spit an error message in case of problem
#define CREATE_COM_INSTANCE(COM_PTR,COM_CLASS) \
{	\
		HRESULT hr;	\
		if (FAILED(hr = COM_PTR.CreateInstance(__uuidof(COM_CLASS)))) {	\
			cerr << "Could not create" #COM_CLASS " instance, hr("	\
							<< hex << hr << dec << ")" << endl;	\
			exit(2);	\
}	\
}

// macro used to setup a pointer
#define CREATE_PTR(COM_PTR,COM_CLASS) \
	I##COM_CLASS##Ptr COM_PTR

// macro used to completely create a C object
// setup a pointer and fill the object
#define CREATE_PTR_AND_COM_INSTANCE(COM_PTR,COM_CLASS) \
	CREATE_PTR(COM_PTR,COM_CLASS);	\
	CREATE_COM_INSTANCE(COM_PTR,COM_CLASS)

#else
/**
comdriver - based on Xiangyang Liu's XYDispDriver

This is a driver that entirely uses COM classes only by using the IDispatch base type.

Thus no need *at all* for any form of automatic COM importation in the compiler (like MSVC's #import),
or even having the COM DLL at compile time to begin with.

Instead, the XYDispDriver searches the type library and calls the methods by name at run time (= late binding).

I have made a "analyst_dispdriver" class&header files which creates Driver classes derivated from the base driver
Each has proper constructor to initialize the IDispatch with the corresponding GUID

 */

// macros used to populate a pointer with an instance of a class
// or spit an error message in case of problem
#define CREATE_COM_INSTANCE(COM_PTR,COM_CLASS) \
	COM_PTR = new COM_CLASS##Driver();

// macro used to setup a pointer
#define CREATE_PTR(COM_PTR,COM_CLASS) \
	COM_CLASS##Driver* COM_PTR

// macro used to completely create a C object
// setup a pointer and fill the object
#define CREATE_PTR_AND_COM_INSTANCE(COM_PTR,COM_CLASS) \
	CREATE_PTR(COM_PTR,COM_CLASS);	\
	CREATE_COM_INSTANCE(COM_PTR,COM_CLASS)

#endif



/*
 * Name spaces
 */

using namespace std;
//using namespace bol;

using namespace ExploreDataObjects;	// object for manipulation of file and data within are here
	// FMANSpecData
	// FMANChromData
	// FMANWiffFile
	// FMANWiffFileControl
using namespace MSMethodSvrLib;	// parameters of experiment are here
	// Experiment

using namespace ExploreDirLib;	// data procession
	// ExploreProcessing 
	// used to perform processing operations (baseline subtract, centroid, noise filter and smooth) directly on a data object.


/*
 * Data types used by Analyst
 */

// scantype
typedef enum scan_enum {
	// from Analyst CookBook
	SCAN_Q1 = 0, QCAN1_MI, SCAN_Q3, SCAN_Q3_MI, SCAN_MRM, SCAN_PREC, SCAN_PROD_ION, SCAN_NEUT_LOSS, 
	// hapilly stolen^H^H^H^H^H^H borrowed from mzWiff
	SCAN_TOF_MS, SCAN_TOF_MS2, SCAN_TOF_PREC, SCAN_EPI, SCAN_ER, SCAN_MS3, SCAN_TDF, SCAN_EMC, SCAN_EMS, SCAN_ELSE, SCAN_UNKNOWN
} SCANENUM;

static const struct {
	const char *name;
	const int level;
} ScanType[] = { 
	// from Analyst CookBook
	{"Q1 Scan",1},
	{"Q1 MI",1},
	{"Q3 Scan",2},
	{"Q3 MI",2},
	{"MRM",2},
	{"Precursor Scan",1},
	{"Product Ion Scan",2},
	{"Neutral Loss Scan",1},
	// shamelessly copied from mzWiff
	{"TOF MS",1},
	{"TOF MS2",2},
	{"TOF Precursor Ion Scan",1},
	{"EPI",2},
	{"ER",1},
	{"MS3",3},
	{"TDF",1},
	{"EMC",1},
	{"EMS",1},
	{"Else",1},
	{"Unknown",1},
	{NULL, 0}};


#define STRINGIFY(x) XSTRINGIFY(x)
#define XSTRINGIFY(x) #x


/*
 NOTE for now the implementation uses flat function and globals
 TODO we could package the globals (with Get/Set -ers), the Analyst object that are passed around, the RunLcmsms (with a Getter) and the functions as an object 
*/


// globals
bool docentroid = false;	// set to true if centroidation must be applied on data
// FIXME program crashes while trying to centroid	

/**
 * dosample processes 1 sample using the set of objects specified by the COM classes
 * @param sample	sample number to process (warning: samples are 1 to N, not 0 to N-1)
 * @param fname	string holding the file name
 * @param wiff	pointer to an instance of coclass FMANWiffFile
 * @param chrom	pointer to an instance of coclass FMANChromData
 * @param spec	pointer to an instance of coclass FMANSpecData
 */
void dosample(const int sample, const char *fname, const char *ofname, 
	CREATE_PTR(wiff,FMANWiffFile),
	CREATE_PTR(chrom,FMANChromData),
	CREATE_PTR(spec,FMANSpecData)) {
/** Hierachical structure of Wiff file :
	- A Wiff file is a collection of data about a series of samples.
	WARNING: The samples are counted from 1 to N (*not* 0 to N-1)

	- A sample can contain several periods - different groups of experiment that were done at different times.
	basically, that means it's possible to do completely different runs on the same sample
	but normally only 1 period is used in most lab workflows.
	WARNING: Whenever looking inside an experiment, Analyst freeze for a certain amount of time processing data
	TODO ask PAB and see with the people of chemical analysis how the elution-time specific detection will be implemented

	- for each period we have several "experiments" : specific settings at which data was collected
	usually most lab workflow use "IDA" :
	- during the whole elution of the sample through the LC column, the spectrometers constantly does a simple analysis of the elute
	(using plain MS or MRM) and continuously store the data.
	this data and the corresponding settings are grouped as "experiment #0"
	- whenever the "experiment 1" picks up something interesting, more detailled analysis will be performed on the elute (some king of MS/MS). 
	Usually these analisys are performed using several different conditions (varying the collision energy).
	each set of data for a given condition is stored in 1 experiment. 
	(most lab workflow will have 3 different collision energies, stored as "experiment #1", "experiment #2" and "experiment #3")

	- the repeated mesures of an experiment are called "cycles". There are as many cycles as there are spectra per experiment or points per chromatogram
*/
	// used to process data - TODO should go as static member once whole reader is objectified
//	CREATE_PTR_AND_COM_INSTANCE(centopts,CentroidOptions);	// TODO change parameters if default aren't ok. Or display Dialog.
//	CREATE_PTR_AND_COM_INSTANCE(process,ExploreProcessing);	// FIXME program crashes while trying to centroid

        RunLcmsms run; // we will store the results in a CaptainBol's runlcmsms structure

        string samname (ConvertBSTRToString(wiff->GetSampleName(sample)));
        cout << "sample " << sample << " - " << samname;

	const int numperiod = wiff->GetActualNumberOfPeriods(sample);
	cout << " - Per:" << numperiod;
	
        run.setDescription(samname);

	for (int pe = 0; pe < numperiod; ++pe) {
		// NOTE Analyst usually freezes here
		const int numexp = wiff->GetNumberOfExperiments(sample, pe); // <- exactly
		// TODO provide a nice status window during this time

		if (numperiod > 1) {
			if (pe == 0) cout << endl;
			cout << "Period:" << pe;
		}
		cout	<< " - Exp:" << numexp;
		
		const int numcycles = wiff->GetActualNumberOfCycles(sample,pe);
		cout	<< " - Cyc:" << numcycles
			<< endl;
		for (int xp = 0; xp < numexp; ++xp) {
		/** as it is an application for controlling an industrial spectrometer doing experiment, interesting data is split across two different place
			- settings, parameters, etc. In short : everything the scientist entered in the software before the machine started to work
			is stored into various structure of the MSMethodSrvLib DLL. 
			for example the collision energy for a given experiment (= series of data, see above) will be stored here.
			- all the actual readings that came from the detection unit of the mass spectrometer, 
			are acessible through the various structure provided by ExploreDataObjects
		*/
		
#ifdef XYDISPDRIVER_H

#if 0
			IDispatch *ptr = static_cast<IDispatch*>(wiff->GetExperimentObject(sample,pe,xp));	// the file it-self
			ExperimentDriver *exp = new ExperimentDriver(ptr);
			
			const SCANENUM scantype = (SCANENUM) exp->GetScanType();
			const int negpolarity = exp->GetPolarity();
			const char *desc = ConvertBSTRToString(exp->GetDescription());
#endif
			// HACK get part of the information from the first spectrum with FMANSpecData->GetScanInfo
			short _gettype = SCAN_UNKNOWN;
			short _getmode = 0;
			long _negative = 0;
			short _dumpq1 = 0;
			short _dumpq3 = 0;
			cout << "T" << endl;
			chrom->SetToTIC(sample,pe,xp);	// look inside the TIC
			float _sec = chrom->GetXValueInSec(1); // get time of fist spectrum
			cout << "S" << _sec << endl;
			spec->SetSpectrum(sample,pe,xp, _sec,_sec); // get the first spectrum
			cout << "I" << endl;
			spec->GetScanInfo(&_gettype, &_getmode, &_negative, &_dumpq1, &_dumpq3);
			const SCANENUM scantype = (SCANENUM)(int)_gettype;
			const int scanmode = _getmode;
			const int negpolarity = _negative;
			// TODO find a way to either use GetExperimentObject to properly initialise the exp object
			const char *desc = ConvertBSTRToString(chrom->GetDataTitle());
#else
			const IExperimentPtr exp = wiff->GetExperimentObject(sample,pe,xp);	// get the settings
			const SCANENUM scantype = (SCANENUM) exp->ScanType;
			const int negpolarity = exp->Polarity;
			const char *desc = exp->Description;
#endif
			// other interesting properties are : ScanMode 

			cout << "\ttype:" << scantype 
				<< "\t" << (negpolarity ? "-" : "+")
				<< " " << ((scantype < SCAN_UNKNOWN) ? ScanType[scantype].name : ScanType[SCAN_UNKNOWN].name) 
				<< "\t\"" << desc
				<< "\"\tlevel:" << ((scantype < SCAN_UNKNOWN) ? ScanType[scantype].level : ScanType[SCAN_UNKNOWN].level);
			cout << endl;

			if (scantype >= SCAN_UNKNOWN) {
				// Unkown : just skip it.
				cout << "\tskipping." << endl;
				continue;
			} 

			/**
			- ChromData object gives a bird's eye view of everything that was hapenning during the whole elution time (chromatogram)
			- SpecData gives detailed data about what was detected at a specific point in time (spectra)
			- At each point in time in a chromatogram, there are spectra aquired
			- we can pick the sum of spectra for a whole time period, but usually we will only want the data at the precise moments of the chromatogram

			WARNING data points are numbered from 1 upward (like sample)
			*/

			chrom->SetToTIC(sample,pe,xp);	// get the data : overall chromatographe
			
			/** IDA
			as explained above LC MS/MS runs are usually setup with a mode that Analyst calls IDA :
			several "experiments" i.e.: several different type of scan are done during the elution.
			
			"experiment #0" - the first type of scan, is systematically done on all cycles
			it is used to quickly identify potential candidate for analysis in the elute.
			often it a simple MS^1 scan - which is simply used to detect precusors and their masses
			it can also be y MRM (see bellow) - with is used to recognise potentially known substance
			
			"experiment #1" and the following - are done only on interesting candidates determined by the previous one.
			they are MS^2 and are used to identify ions found on an MS^1
			or to analyse and confirm identity of potential candidates found with MS^2
			*/
			
			if (scantype == SCAN_MRM) {
				/** MRM - Multiple reaction monitoring

				mostly like Reaction Monitoring (a single peak MS^2) which is made to quickly identify candidates that are potentially known substances :
				- based on their mass (q1 with fixed mass)
				- and on a specific fragment (q3 with fixed mass)
				
				an MRM is a collection of such Reaction monitoring, designed to look for multiple such candidates
				- it could be considered as a collection of multiple single-peak MS^2
					(that would be the most logical for bioinformatics processing)
				- or as a weird spectrum whose X axis aren't simply M/Z, but M/Z pairs (for q1 and q3)
					(that's what Analyst does. It has the advantage of being more compact to display on screen)
				
				This type of scan is useful when screening a sample for known substances.
				
				HACK currently we consider that an MRM shows precursors
				so the MRM is processed to create a "virtual" MS^1 spectrum of precursors
				it's not fundamentally wrong, as traditionnally those are represented as chromatograms
				(per each precursors, a collection of intensity over time)
				wheras we simply have the data packaged the other around
				(per each time slot, a collection of intensity over precursors' m/z)

				TODO at some time later once the mecanism is correctly added in CaptainBOL, we should use a special calsse
				PeakIntensityRM which can also hold additionnal informations about the MRM/SRM
				*/
				cout << "\tprocessing MRM..." << endl;
				if (xp > 0) 	cout << "\tshould be the first ?!?";
				cout << "\tgetting precursor..." << endl;
				
				const int chpts = chrom->GetNumberOfDataPoints();
				cout << "\tNumber of data points: " << chpts;

				double min, max;
				chrom->GetXValueRange(&min, &max);
				cout << "\tThere's data between " << min << " and " << max << endl;

				// scan the whole chromatogram -> the whole duration of the elution
				for (int j = 1; j <= chpts; ++j) {	// Warning: DataPoints are 1 to N (not 0 to N-1 like cycles)
				// this loop parse the mass readings
				/** semantic consideration about MS measurment
				
				- inside a chromatogram, it's simply called a "data point" - Beware: 1 to N (*not* 0 to N-1) this contradicts "cycles"
				- inside the "msmethod" classes, it's called a cycle - Beware : 0 to N-1 (*not* 1 to N)
				- spectra aren't addressed directly but by time-frame - Beware : this is *always* specified in seconds (whatever the default time scale of the chromatogram is)
				*/
					const double z = chrom->GetDataPointYValue(j); // total intensity at this point in time / this cycle
	
					// are there at least some detected ions ?
					if (! (z > 0.)) continue; 
					
					// (So we'll use the time reported by the chromatogram to get corresponding spectra)
					const double t = chrom->GetXValueInSec(j); // we use seconds as standard [even if some wiff file are produced using minutes]

					cout << "\t\tProcessing MRM " << j << "\tat " << t;

					spec->SetSpectrum(sample,pe,xp,static_cast<float>(t),static_cast<float>(t)); // at the same exact time as the point on the chromatogram
						// note : whatever the default time scale in the wiff file, the "SetSpectrum" method *needs* seconds.

					const int sppts = spec->GetNumberOfDataPoints();

					const char *cptr = ConvertBSTRToString(spec->GetDataTitle());
					const string title(cptr);
					delete cptr;

					cout	<< "\tnum of points: " << sppts
						<< "\tspectra:" << spec->GetNumberOfScans()
						<< endl;
					// NOTE no centroiding as MRM aren't continuous
	
					/*
					 * HACK sum together intensities for precursors with similar mass
					 * TODO once PeakIntensityRM are implemented this should be deprecated and we simply build an instance of PeakIntensityRM using Q1 and Q3
					 */

/* TODO once PeakIntensityRM are implemented 
					SpectrumLcms* spectrum=new SpectrumLcms(t);
					spectrum->setDescription(title);
					for (int i = 1; i <= sppts; ++i) {
						double pit = spec->GetDataPointYValue(i);
						
						if (! (pit > 0.)) continue; // skip if there's no peak intensity
						double pmz = 0., fmz = 0.;
						spec->GetQ1Mass(i, &pmz);
						spec->GetQ3Mass(i, &fmz);

						spectrum->push_back(new PeakIntensityRM(pmz, fmz, 0, pit));
					}

					// congratulation, we'd created a virtual MS^1 spectrum
					run.msCollection().push_back(spectrum);
*/
					// first cache all Q1 (=precusor) mass settings
					// the Q3 (=fragment) are simply ignored
					double *mz = new double [sppts];
					double *it = new double [sppts];
					int csize = 0;
					for (int i = 1; i <= sppts; ++i) {
						double pit = spec->GetDataPointYValue(i);
						
						if (! (pit > 0.)) continue; // skip if there's no peak intensity
						
						double pmz = 0.;
						spec->GetQ1Mass(i, &pmz);
						mz[csize] = pmz;
						it[csize] = pit;
						++csize;
					}

					SpectrumLcms* spectrum=new SpectrumLcms(t);
					spectrum->setDescription(title);
					for (int i = 0; i < csize; ++i) {
						// skip flagged M/Z
						if (mz[i] <= 0.) continue;
						
						double	pmz = mz[i],
      							pit = it[i];
						// sum all mass close to pmz
						for (int k=i+1; k < sppts; ++k) {
							if (fabs(pmz - mz[k]) < 0.01) {	// ridiculously close delta. TODO maybe DeltaMass should work better here...
								mz[k] = 0.;
								pit += it[k];
							}
						}
						
						// add the virtual MS^1 peak to the spectum
						if (pit > 0.) // (only if there's peak intensity)
							spectrum->push_back(new PeakIntensity(pmz, 0, pit));
					}
	
					// congratulation, we'd created a virtual MS^1 spectrum
					run.msCollection().push_back(spectrum);
					
					delete mz; delete it;
				} // cycles / points on the chromatogram

				continue;
			} else if (ScanType[scantype].level == 1) {
				// MS level 1 - usually the experiment #0 in an IDA
				// used to detect precursors
				if (xp > 0) 	cout << "\tshould be the first ?!?";
				cout << "\tgetting precursor..." << endl;
				
				const int chpts = chrom->GetNumberOfDataPoints();
				cout << "\tNumber of data points: " << chpts;

				double min, max;
				chrom->GetXValueRange(&min, &max);
				cout << "\tThere's data between " << min << " and " << max << endl;

				// scan the whole chromatogram -> the whole duration of the elution
				for (int j = 1; j <= chpts; ++j) {	// Warning: DataPoints are 1 to N (not 0 to N-1 like cycles)
				// this loop parse the mass readings
				/** semantic consideration about MS measurment
				
				- inside a chromatogram, it's simply called a "data point" - Beware: 1 to N (*not* 0 to N-1) this contradicts "cycles"
				- inside the "msmethod" classes, it's called a cycle - Beware : 0 to N-1 (*not* 1 to N)
				- spectra aren't addressed directly but by time-frame - Beware : this is *always* specified in seconds (whatever the default time scale of the chromatogram is)
				*/
					const double z = chrom->GetDataPointYValue(j); // total intensity at this point in time / this cycle
	
					// are there at least some detected ions ?
					if (! (z > 0.)) continue; 
					
					// (So we'll use the time reported by the chromatogram to get corresponding spectra)
					const double t = chrom->GetXValueInSec(j); // we use seconds as standard [even if some wiff file are produced using minutes]

					cout << "\t\tGetting MS1 spectrum " << j << "\tat " << t;

					spec->SetSpectrum(sample,pe,xp,static_cast<float>(t),static_cast<float>(t)); // at the same exact time as the point on the chromatogram
						// note : whatever the default time scale in the wiff file, the "SetSpectrum" method *needs* seconds.
					
					bool centroided = spec->GetIsCentroided(); 
					
					if ((! centroided) && docentroid) {
//						process->Centroid(centopts,spec);
						centroided = true;
					}						

					const int sppts = spec->GetNumberOfDataPoints();

					const char *cptr = ConvertBSTRToString(spec->GetDataTitle());
					const string title(cptr);
					delete cptr;
#ifdef XYDISPDRIVER_H
					int stype = spec->GetSpectrumType();
#else
					SpectrumType stype;
					spec->get_SpectrumType(&stype);
#endif


					cout	<< "\tnum of points: " << sppts
						<< "\tspectra:" << spec->GetNumberOfScans()
						<< "\t" << (centroided  ? "centroided" : "raw")
						<< endl;

	
					SpectrumLcms* spectrum=new SpectrumLcms(t);
					spectrum->setDescription(title);

					for (int i = 1; i <= sppts; ++i) {
						double pmz, pit;
						spec->GetDataPoint(i, &pmz, &pit);
				
						if (pit > 0.) // (only if there's peak intensity)
							spectrum->push_back(new PeakIntensity(pmz, 0, pit));
					}

					/* // currently ranks aren't needed in a loader

					SpectrumPeakIntensityUtilities::computeIntensityRank(*spectrum);
					sort(spectrum->begin(), spectrum->end(), Peak::ltptr);

					// to get spectra numbers
					cout << run.msCollection().size() << endl;
					*/
					run.msCollection().push_back(spectrum);
				} // cycles / points on the chromatogram

				continue;
			} else if (ScanType[scantype].level == 2) {
				// experiment #1 and following : detailed analysis
				if (xp == 0) cout << "\tshould be level 1 ?!?";
				if (run.msCollection().empty()) {
					cerr << "Ooops: no MS^1 from which to pick percursors" << endl;
				}
				
				
				cout << "\tgetting spectra..." << endl;

				const int chpts = chrom->GetNumberOfDataPoints();
				cout << "\tNumber of data points: " << chpts;

				double min, max;
				chrom->GetXValueRange(&min, &max);
				cout << "\tThere's data between " << min << " and " << max << endl;

				// scan the whole chromatogram -> the whole duration of the elution
				for (int j = 1; j <= chpts; ++j) {	// Warning: DataPoints are 1 to N (not 0 to N-1 like cycles)
				// this loop parse the mass readings
				/** semantic consideration about MS measurment
				
				- inside a chromatogram, it's simply called a "data point" - Beware: 1 to N (*not* 0 to N-1) this contradicts "cycles"
				- inside the "msmethod" classes, it's called a cycle - Beware : 0 to N-1 (*not* 1 to N)
				- spectra aren't addressed directly but by time-frame - Beware : this is *always* specified in seconds (whatever the default time scale of the chromatogram is)
				*/
					const double z = chrom->GetDataPointYValue(j); // total intensity at this point in time / this cycle
	
					// are there at least some detected ions ?
					if (! (z > 0.)) continue; 
					
					// (So we'll use the time reported by the chromatogram to get corresponding spectra)
					const double t = chrom->GetXValueInSec(j); // we use seconds as standard [even if some wiff file are produced using minutes]

					cout << "\t\tGetting spectrum " << j << "\tat " << t;

					spec->SetSpectrum(sample,pe,xp,static_cast<float>(t),static_cast<float>(t)); // at the same exact time as the point on the chromatogram
						// note : whatever the default time scale in the wiff file, the "SetSpectrum" method *needs* seconds.
					
					bool centroided = spec->GetIsCentroided();
					if ((! centroided) && docentroid) {
//						process->Centroid(centopts,spec);
						centroided = true;
					}
					

					const int sppts = spec->GetNumberOfDataPoints();
					double it = z, // initialise intensity with some bogus data
						MoZ = 0.0, m2 = 0.0, ce = 0.0; // ce = collision enery
					short type = SCAN_UNKNOWN, ch = 0;	// ch = charge state

					// NOTE the following Method isn't documented in the CookBook, but was obtained by speaking with the technical support of Applied
					wiff->GetScanInfo(sample,pe,xp,j-1,// Warning: cycles are 0 to N-1 (not 1 to N like DataPoints)
						&type,&MoZ,&m2,&ce,&ch);

					cout << "\tce:" << ce << " ch:" << ch  << " t:" << type;

					// find which was the precursor's intensity in the MS^1 collection
					//	- at the SAME OR PREVIOUS time
					//	- at m/z "MoZ" within a small delta
					SpectrumPeakRefDualPtr<PeakIntensity>* refptr = NULL; // this will hold the reference from the MS^2 spectrum to the MS^1 and peak
					if (! run.msCollection().empty()) {
						it = 0;
						// scan time axis
						// TODO this will be slow on big runs. A timeindex (similar to the mass index for MS^2) could be handy
						
						// NOTE this works because analyst stores the spectra ordered by time inside a given experiment
						// and thus, the spectra where added to RunLcmsms chronologically.
						
						// first find the last MS^1 before current t
						SpectrumLcmsCollection::iterator time = run.msCollection().begin();
						// post-fix : the *previous* loop was the last. (i.e.: do one more loop *with* the tail)
						while ((time++) != run.msCollection().end()) 
							if ((time == run.msCollection().end()) // past the last ? the last is the one we need.
							   || (*time)->retentionTime() > t) { // found 1 which is older than t
								if (time == run.msCollection().begin())
									break; // even the first MS^1 is older than t ? Bail out !
								else
									--time; // move back to the last whose retention time was befor curent t.
								
								// Ok, we find the preceding MS^1 spectum
								// now scan all masses
								SpectrumLcms::iterator mass = (*time)->begin();
								while (mass != (*time)->end()) 
									if (fabs(MoZ - (*mass)->moz()) < .1) { // TODO use a DeltaMass
										// Match ! store the intensity
										it = (*mass)->intensity();
										// and build the reference
										refptr =  new SpectrumPeakRefDualPtr<PeakIntensity> ((*time), (*mass));
										break;
									} else
										++mass;
								break;
							}
					} else {
						it = z; // no MS^1, put some bogus value based on total ions
					}

					const char *cptr = ConvertBSTRToString(spec->GetDataTitle());
					const string title(cptr);
					delete cptr;
					
#ifdef XYDISPDRIVER_H
					int stype = spec->GetSpectrumType();
#else
					SpectrumType stype;
					spec->get_SpectrumType(&stype);
#endif


					cout << "\tPre: " << MoZ << "\tInt: " << it 
						<< "\tnum of points: " << sppts
						<< "\tspectra:" << spec->GetNumberOfScans()
						<< "\ttype:" << stype
						<< "\t" << (centroided  ? "centroided" : "raw")
						<< endl;

					// TODO ask PAB what are the usual range of charges and see if "charge state" above is revelent
	
					SpectrumLcmsms* spectrum=new SpectrumLcmsms(MoZ, Peak::string2chargemask("2+ or 3+"), it);
					spectrum->setDescription(title);
					spectrum->setRetentionTime(t);
					spectrum->setCollisionEnergy(ce);
					if (refptr != NULL) {
						// we have a reference to a MS^1 spectrum & peak....
						SpectrumPeakRefPtr baseptr(refptr);
						spectrum->setSpectrumPeakRef(baseptr);

						// if everything is okay, now we can travel from the current spectrum to the parent intensity
						// (and both should be the same)
						assert(spectrum->spectrumPeakRef()->peak()->intensity() == it);
					}
					
					for (int i = 1; i <= sppts; ++i) {
						double pmz, pit;
						spec->GetDataPoint(i, &pmz, &pit);
				
						if (pit > 0.) // (only if there's peak intensity)
							spectrum->push_back(new PeakIntensity(pmz, 0, pit));
					}

					/* // currently ranks aren't needed in a loader

					SpectrumPeakIntensityUtilities::computeIntensityRank(*spectrum);
					sort(spectrum->begin(), spectrum->end(), Peak::ltptr);

					// to get spectra numbers
					cout << run.msmsCollection().size() << endl;
					*/
					run.msmsCollection().push_back(spectrum);
				} // cycles / points on the chromatogram

				continue;
			} else if (ScanType[scantype].level == 3) {
				// experiment #1 and following : detailed analysis
				if (xp == 0) cout << "\tshould be level 1 ?!?";
				cout << "\tgetting MS3..." << endl;	
				
				// TODO implement MS^3 in RunLcmsms
				cout << "\t\tNot implemented yet !!!" << endl;	
				continue;
			}

		} // experiment

	} // period
	cout << endl;
	cout << endl;
	
	boost::regex rx_remove_extension("\\.[^./\\\\]*$", boost::regex_constants::perl);
	boost::regex rx_remove_nonalphanum("[^-_.~[:alnum:]]", boost::regex_constants::perl);

	string inname (fname);
	string samsuf (boost::regex_replace(samname,rx_remove_nonalphanum,"", boost::match_default | boost::format_all));
#if 0
	{
          string newsuf ("-" +  boost::lexical_cast<string>(sample) + "-boostser.xml");

          string outname (boost::regex_replace(inname,rx_remove_extension,newsuf));
          const char *cptr = outname.c_str();
          cout << "outname:\t" << cptr << endl;

#ifdef HAVE_BOOST_SERIALIZATION
          // save data to archive

          ofstream ofs(cptr);
          boost::archive::xml_oarchive oa(ofs);
          oa << BOOST_SERIALIZATION_NVP(run);
#else
          cerr << "boost serialisation isn't compiled in" << endl;
#endif
        }
#endif

	{
		string outname (ofname
			? ofname 
			: (boost::regex_replace(inname,rx_remove_extension,"") + "-" +  boost::lexical_cast<string>(sample) + "-" + samsuf + ".sml"));
	
		const char *cptr = outname.c_str();
		cout << "outname:\t" << cptr << " ..." << endl;

#ifdef HAVE_PROTOBUF
		// save data to archive
		fstream os(cptr, ios::out | ios::trunc | ios::binary);
	
		protobuf::RunLcmsms pb;
	
		// operator declared in the include file will take care of populating the google protobuf serializer
		pb << run;
	
		if (!pb.SerializeToOstream(&os)) {
			cerr << "Failed to serialize data." << endl;
		}
		cout << "done." << endl;
#else
		cerr << "google protocole buffer serialisation isn't compiled in" << endl;
#endif
	}
	// to manually inspect output :
	//  protoc --decode=bol.protobuf.runlcmsms runlcmsms.proto < /tmp/testser.bin  | less

        return;
}



/**
 * dofile takes a Wiff file and process its sample
 * @param _fname : path to a WIFF file
 * @param sample : sample number to process. By default : 0 = process all samples from file
 */
void dofile(const char *_fname, const char *ofname, int sample = 0) {
	BSTR fname = ConvertStringToBSTR(_fname);

	cout << "reading from : " << fname << endl;
/*
	NOTE Normally, according to the Analyst Cookbook, we open a file using the FileControl with something like this
	CREATE_PTR_AND_COM_INSTANCE(filecontrol,FMANWiffFileControl);
	IFMANWiffFilePtr wiffile = filecontrol->GetWiffFileObject(fname, 1)

	*BUT* surprisingly, lots examples in theCookBook use directly FMAN*Data objects,
	using  "->WiffFileName" property to get/set the filename
	and "->GetWiffFileObject" to get the wiff file object.

	mzWiff uses the same approach
*/

	CREATE_PTR_AND_COM_INSTANCE(chrom,FMANChromData);
	CREATE_PTR_AND_COM_INSTANCE(spec,FMANSpecData);

	// set-up objects for interacting with file
#ifdef XYDISPDRIVER_H
	chrom->PutWiffFileName(_fname);	// for chromatogram
	spec->PutWiffFileName(_fname);	// for specta
	

	IDispatch *ptr = static_cast<IDispatch*>(chrom->GetWiffFileObject());	// the file it-self
	cout << ptr << endl;
	// For some reason, Attaching to the output of GetWiffFileObject fails with XYDispDriver
	cout << "Attach... "  << endl;
	FMANWiffFileDriver *wiff = new FMANWiffFileDriver (ptr);
	cout << "Attached"  << endl;
	
/*
	// the FMANWiffFileControl isn't accesible either.with XYDispDriver
        // maybe see this ?
        XYDispDriver disp;
        if(disp.CreateObject(_T("Word.Application"), CLSCTX_LOCAL_SERVER))

        By default, the second parameter is CLSCTX_ALL, it doesn't work for some types of objects on some platforms.
*/
	
/*
	// HACK we will use the OpenWiffFile method instead
	// TODO find a way to either use FMANWiffFileControl or GetWiffFileObject to properly initialise the wiff object
	CREATE_PTR_AND_COM_INSTANCE(wiff,FMANWiffFile);
	wiff->OpenWiffFile(_fname);
*/
#else
	chrom->PutWiffFileName(fname);	// for chromatogram
	spec->PutWiffFileName(fname);	// for specta
	IFMANWiffFilePtr wiff = chrom->GetWiffFileObject();	// the file it-self
#endif

	int numsamples = wiff->GetActualNumberOfSamples();
	cout << "num of sample :\t" << numsamples << endl << endl;
	
	if (sample > 0) {
		dosample(sample	, _fname, ofname, wiff, chrom, spec);	
	}else {
		// process all samples in the file
		for (int i = 1; i <= numsamples; i++)	// Warning: Samples are 1 to N (and not 0 to N-1 like mostly everything else)
			dosample(i, _fname, ofname, wiff, chrom, spec);
	}

	SysFreeString(fname);
}





int main(int argc, char *argv[]) {
	cout << "CaptainBol is in da house !!!" << endl << endl;
	
	using namespace boost::program_options;
	int sample = 0;
	string ifname ("");
	string ofname ("");
	
	// Declare the supported options.
	options_description desc("Allowed options");
	desc.add_options()
		("help,h", "produce help message")
		("sample,s", value<int>(&sample)->default_value(0), "proceed one single sample")
		("input-file,i", value< string >(&ifname), "input file name")
		("output-file,o", value< string >(&ifname), "output file name for one single sample")
		;
	
	positional_options_description pos;
	pos.add("input-file", 1);
	pos.add("output-file", 1);

	variables_map vm;

	store(command_line_parser(argc, argv).options(desc).positional(pos).run(), vm);
	notify(vm);

	if (! vm.count("input-file")) {
		cerr << "Please at least give an input file name" << endl
				<< desc << endl;
		exit(1);
	}
	if (vm.count("help")) {
		cout << desc << endl;
		return 1;
	}

	if (sample) {
		cout << "Processing sample " 
			<< sample <<  endl;
	} else {
		cout << "Process all samples" << endl;
		if (vm.count("output-file")) {
			cerr << "One output file can only hold one single sample" << endl
				<< "- either specify sample number to process" << endl
				<< "- or don't specify an output file name and let the software generate names" << endl
					<< desc << endl;
			exit(1);
		}
	}
	

#ifndef XYDISPDRIVER_AUTO_INIT
/**
 * Never forget to initialise first otherwise "800401f0 - CoInitialize has not been called."
 */

//	switch (CoInitializeEx(NULL, COINIT_SPEED_OVER_MEMORY|COINIT_APARTMENTTHREADED)) { // COINIT_APARTMENTTHREADED COINIT_MULTITHREADED
	switch (CoInitialize(NULL)) {
		case S_OK:
			cout << "COM initialised" << endl;
			break;
		case S_FALSE:
			cout << "COM already initialised in this thread" << endl;
			break;
#ifdef RPC_E_CHANGED_MOD
		case RPC_E_CHANGED_MOD:
			cerr << "different concurrency model  or neutral threaded apartment" << endl;
			exit(2);
#endif
		case E_INVALIDARG:
			cerr << "invalid arg" << endl;
			exit(2);
		case E_OUTOFMEMORY:
			cerr << "out of memory" << endl;
			exit(2);
		case E_UNEXPECTED:
			cerr << "unexpected error" << endl;
			exit(2);
		default:
			cerr << "unkown COM initialisation failure" << endl;
			exit(2);
	};
#endif

	dofile(ifname.c_str(), vm.count("output-file") ? ofname.c_str() : NULL, sample);

	// all instances of COM objects are handled with smart pointer and should automatically collect garbage
	CoUninitialize();
	return EXIT_SUCCESS;
}


// CODE GRAVEYARD :

/* // not used anymore
// Until we find a better way, this regular expression extracts the precursor Mass from the title of the spectrum
boost::regex re_find_prec("^[^\\(]+\\(([^\\)]+)\\).*$", boost::regex_constants::perl); 

						// no information in documentation about how to get a spectrum's precursor. For now, we are parsing the title. Yes, indeed. Kill me for it.
						boost::cmatch matches;
						boost::regex_match(cptr, matches, re_find_prec);
						string found(matches[1]);
						MoZ = boost::lexical_cast<double>(found);
						// hope that one day Applied will answer my support mail and tell how I'm actually supposed to do...
*/
